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Abstract 

We propose a new phenomenology of the generation of internal 
transport barriers, based on the exact periodic solution of the Flierl- 
Petvishvili equation. We examine the stability of this solution and 
compare the late stages of the flow with the ensemble of vortices. 
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1 Introduction 

The aim of these notes is to contribute to the clarification of the dynamical 
processes leading to the formation, confinement character sit ics and intermit- 
tent behavior of Internal Transport Barriers (ITB). We focus here on the 
role of the stationary periodic structure of the flow appearing in the inter- 
mediate spatial range between the highly radially elongated eddies of the Ion 
Temperature Gradient (ITG) instability and the Larmor - radius scale vor- 
tices described by the Hasegawa-Mima equation. In this "mesoscopic" range 
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the structures are described by the ion dynamics equation where the scalar 
nonhnearity (in contrast to the polarisation drift nonlinearity) is dominating. 

The current physical understanding of the formation of the transport 
barriers is based on the effect of the sheared poloidal plasma rotation on 
the ITG potential structures. The ITG potential structures are in general of 
large radial extension. They may persist even in the turbulent regime and 
are a major agent of transport since the plasma convection inside the eddies 
is an efficient support for temperature advection between high (small r) and 
low (high r) temperature regions. It is well known that this geometry of the 
flow generates velocity stress and that the divergence of this tensor induces a 
substantial contribution in the momentum balance equation for ions. When 
this momentum drive is sufficiently high, the plasma begins to rotate, and 
a space distribution of the poloidal velocity is established. The shear of the 
poloidal velocity has a strong effect on the turbulent structures with com- 
parable space extension and the radial correlation length is reduced. This 
correspondingly reduces the transport leading to the formation of a transport 
barrier. A theory that explains this in quantitative terms has been proposed 
by Diamond and coworkers, as a "predator-pray" model, in general formal- 
ized with Fischer-type equations. It consists of the description, in terms of 
the density of the energy quanta, of the balance in phase space of two distinct 
populations, associated with the turbulent waves and respectively with the 
zonal flows. The model naturally leads to the self-organised state consisting 
of mutual control of these two physical entities, found in a permanent rela- 
tive adjustement. This model is essentially a thermodynamical model, which 
probably remains true for a wide variety of dynamical details. 

A picture of the dynamics is allways useful, even if the scaling laws (of 
threshold and transport rate) are in general based on thermodynamical bal- 
ance and conservation properties. The later is a more general (and safer) 
level of description, leading to a picture where the invariance properties of 
the dynamical equations govern the scaling. The dynamical evolution of the 
system is a detailed description and its validity is strongly dependent on 
a clear separation of space-time scales and a correct choice of initial and 
boundary conditions. However, the physical content can only revealed by 
the detailed description. 

In this work wc develop a possible phenomenological model of generation 
of the Internal Transport Barriers in tokamak. This model relies on the role 
of the stationary periodic exact solution of the Flierl-Petviashvih equation as 
an attracting state for the ion turbulence, when the regularity of the eddies is 
still pronounced. We also analyse analytically and numerically the stability 
of this solution. Some properties are shown to be connected with the fact 
that this equation is closely related with other, exactly integrable models. 
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2 Possible phenomenology of the ITB forma- 
tion 



Numerical simulations show that the structure of the potential perturbation 
in the ion temperature gradient (ITG) instability are elongated and rather 
thin on the direction that is transversal to their longest dimension. The 
Reynolds stress arising from the broken symmetry configuration induces the 
rotation of the plasma. It may be useful here to recall the similar situation 
arising in the thermal convection associated with the Rayleigh-Benard (RB) 
instability. Experiments carried out by Howard and Krishnamurti in typical 
RB geometry at increasing Rayleigh number have revealed that the classical 
conduction-convection bifurcation is only the first of a series of transitions 
generating new system 's behaviors. The RB geometry of convective rolls is 
sensitive to perturbations and this leads to a symmetry breaking consisting 
of deformations of the rolls in direction parallel with the two plates, in one 
direction for the low (high temperature) plate and opposite direction for 
the upper (low temperature) plate. This is accompanied by a nonvanishing 
statistical average of the velocity stress tensor, whose nonzero divergence 
(Reynolds stress) is equivalent to a momentum drive. The lower part of 
the convective fluid is entrained in a flow in one direction (parallel to the 
plates) while the upper part of the convective field aquires a flow motion 
in the opposite direction. These flows are called "winds" and they are seen, 
together with the strong deformation of the rolls, on the experimental picture 
of the fluid section. In addition, intermittent processes consisting of rising of 
flnite volumes of cold fluid (from the bottom region) toward the top plate and 
fall of warm fluid from the top region into the cold region have been observed. 
These events are called "plumes" . One may consider from this picture only 
the part relating to the effect of the Reynolds stress as particularly relevant 
for the plasma convection in the eddies of the ITG instability. 

The result of the strong Reynolds stress arising in the ITG instability is 
the deformation ("tilting") of the convection eddies and the generation of 
the corresponding winds. The convective structures are strongly elongated 
and tend to align with the poloidal direction, as imposed by the momentum 
from the Reynolds stress. At a certain moment the deformed ITG pattern 
of flow evolves to the solution consisting simply of parallel layers of flow, a 
structure which is essentially poloidally oriented and periodic along the radial 
direction. We have proved that this is an exact and robust solution of the ion 
dynamics on mesoscopic space scales, (approx. Ps/s) where it is dominated 
by the scalar nonlinearity. Attaining this state, the plasma has practically 
vanishing transport in the radial direction. We have to note that it is not 
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necessary to invoke the tearing appart of the ITG eddies and the destruction 
of the radial correlation via sheared flow. The solution intrinsically has a 
radial scale and no transport because there is no phase mismatch between 
particles and potential. 

The flow represented by this solution may last indefinitely except for the 
coUisional decay which however may be considered low in high temperature 
regimes. Furthermore, in the ideal stationary case, the Reynolds stress van- 
ishes due to the poloidal symmetry. However there is a hmitation of the 
persistance of this periodically layered flow and this arises from the stabil- 
ity against perturbations. We can suppose that this may take place in two 
different ways. The first of them is essentially a reversed process as that 
which led to the periodic layers, and consists of closing the flow lines at finite 
poloidal wavelength, with suppression of the zonal flow and reformation of 
the tilted ITG pattern. The finite radial projection leads to a siidden increase 
of transport. The other way is a strong qualitative change of the flow. 

We have found that the stability of the periodic flow pattern is deter- 
mined both by the amplitude and by the geometry of the initial perturbation. 
Some perturbations (like, for example, the monopolar vortices embedded in- 
side a layer) have a very long stability time since they accomodate with the 
background flow by reshaping the distribution of local velocities. However, 
perturbations that have a high amplitude relative to the background flow 
and/or do not conform to the flow geometry lead to the destabilisation and 
eventually destruction of the flow pattern. It is important however to note 
that the numerical simulations (presently with limitted precision) and ana- 
lytical considerations suggest that the destruction of the flow pattern does 
not immediately lead to an arbitrary random field: a fundamental process is 
the generation of small (~ Ps) space scale vortices. We can show that the 
regular flow structure is replaced by a lattice of vortices which is reminescent 
of the exact solution of a closely related nonlinear equation. This solution 
evolves, due to the weak interaction between the vortices, to an ensemble of 
quasi-independent vortices that collide inelastically and become of various 
amplitudes plus a surrounding drift wave radiation. Prom this random field 
the ITG instability may regenerate the structure of eddies and reinstate the 
transport. The time scale for doing this is ^yiTG, followed by saturation, as 
shown by many simulations. 

The successive steps of this qualitative scenario can repeate themselves 
over and over but the senstivity to perturbation of the time scales involved 
may lead to an intermittent (bursty) behavior rather than to a limit cicle. 
Prom purely theoretical point of view it is hard to expect that this evolution 
can be lumped into a single equation whose phase space could provide the 
successive transitions. 
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We note that the transport is low in the phase where the flow is described 
by the stationary periodic solution of the Flierl-Petviashvili equation and in 
the state where the ensemble of vortices has replaced it (since the elementary 
space scale is very small, of the order of few Larmor radius, as we will show 
below. The large eddies have a substantial transport. 



3 Derivation of the equation at the interme- 
diate space scales 

3.1 The equation for the two-dimensional ion dynam- 
ics 

Consider the equations for the ITG model in two-dimensions with adiabatic 
electrons: 

^ + V-(vinO = 

, e , „ e „ 

^+ v,-V = — -V0 + — ViXB 
at rrii rrii 

We assume the quasineutrality 

rii ^ n,, 

and the Boltzmann distribution of the electrons along the magnetic field line 

He = no exp I — — 

In general the electron temperature can be a function of the radial variable 

The velocity of the ion fluid is perpendicular on the magnetic field and is 
composed of the diamagnetic, electric and polarization drift terms 

T, 1 dn,^ 



e\ B Hi dr ^ 
— V0 X n 



B 

1 f d 



BVti \dt 
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Introducing this velocity into the continuity equation, one obtains an equa- 
tion for the electrostatic potential (f). 

In the Appendix it is presented a derivation of three versions of the 
nonlinear differential equations that may govern, in certain conditions, the 
two-dimensional ion dynamics at intermediate scales. 



3.2 The stationary equation at intermediate space scales 
(Flierl-Petviashvili equation) 

The ion drift wave equation has a distinct dynamical behavior according 
to the space-time scales involved. The Hasegawa-Mima-Charney (HMC) 
equation is obtained for small scales and the stationary states exhibit dipolar 
structures of the order of the Larmor radius. On larger scales (intermediate 
between the HMC dipole and the large ITG eddies) the scalar (or KdV-type) 
nonlinearity is prevailing and the only known structures are monopolar jlj. 
Both are not solitonic but very robust and long lived. In the latter case a one- 
dimensional version of the equation has been derived by Petviashvili |Sj who 
also used it in the study of the Jupiter's Red Spot. The equation has been 
rederived along with a careful analysis of the scales involved IE] , [El , resolving 
a controversy on the role of the temperature gradient [Oj, [ZI. In the study of 
ocean flows, Flierl ITU] has independently formulated an equation with the 
same structure. The one dimensional equation has been solved on an infinite 
domain ^1], [Hj, obtaining as solution the KdV soliton. The two dimensional 
equation has vortical monopolar solutions, well studied numerically [12] , [E] . 
Other applications include long waves on thin liquid films and Rossby waves 
in rotating atmosphere. 

The Flierl-Petviashvili equation (derived several times in the Appendix) 

is 

A0 = a0 - 



where 



Te d fl_\ ^ e d f]_ 
2u'^eBlp'ldx \LnJ 2miu'^ dx \L„ 

By defining a functional of the solution expressed in one-dimensional (ra- 
dial) geometry and taking the extremum of the functional under the condition 
ofasymptotic decay, Petviashvili and Pokhotelov [T^ have found the solution 



4.8a 

mir) = 

' 2(3 



sech ( - yx^ + (y — ut) r 



4/3 

2 
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This is only an approximate solution and the reason to look for such form 
resides in the previous indication that the one dimensional version of this 
equation, with generalized nonlinearity, presents exact solutions as powers of 
the sech function. 

3.3 The monopolar vortex and the similarity with in- 
tegrable equations 

It is well known that the Flierl-Petviashvili equation has long been stud- 
ied in relation with the possibility that it has as solution an isolated vortex 
with nontrivial stability properties. This vortex has never been determined 
analytically and the explanation seems to be that the equation is not inte- 
grable. However, analytical approximations exists Eq.Q and also a very well 
documented study of numerically-calculated vortex solution is available [T7] , 

nni. 

The monopolar vortex has been considered in relation with the Great Red 
Spot, a vortical structure with a long time life (more than three hundread 
years) , which is observed to be embedded in the zonal flow of the atmosphere 
of Jupiter. The condition imposed to the Flierl-Petviashvili equation is to 
provide an isolated, finitely extended vortical solution, obeying boundary 
conditions at infinity (on the 2D) of smooth decaying. 

We consider that the insuccess in determining the analytical vortex solu- 
tion has a fundamental motivation. 

First of all the equation does not pass the Painleve test, so it is very 
likely that it is not integrable. This, in the definition that is suggested by 
the Inverse Scattering Transform method, means that there is no Lax pair of 
operators for this equation. In other words, one cannot find a system of two 
linear differential equations, whose condition of compatibility to be expressed 
as the nonlinear FP equation. 

This equation has close ressemblance with some approximative form of 
other nonlinear equations. For example, consider the equation 

A(f) = exp (—0) 

This is the Liouville equation, with numerous applications. In particular 
it has been invoked in the description of the distribution of current in the 
magnetic confinement systems (tokamak). This equation is integrable and 
exact solutions are known. In addition, exact analytical solutions can be 
constructed on periodic domains. The expansion of the right hand side of 
this equation takes a form close to the FP equation. 



8 




Figure 1: The exact solution of the Jacobs- Rebbi (Abehan-Higgs) equation 
for ideal fluids obtained analytically as a ratio of Riemann B functions. 



Even more close is the differential equation 

A0 = exp (0) - 1 

sometimes known as Abelian-Higgs equation (AH), governing the vortices of 
the superconducting media. This equation has also been derived by Jacobs- 
Rebbi. It is clear that an expansion of the right hand side leads to something 
close to the FP equation, especially because the FP equation has actually 
been derived under the neglect of the third order powers of 0. 

= exp (0) - 1 




Or, this equation is exactly integrable on periodic domains. We have proved 
jTHj that it possesses a Lax pair of operators and we have constructed exact, 
analytical solutions, expressed as ratios of Riemann G functions. 

The solutions of the AH equation, as shown in the Fig^have the config- 
uration of a lattice of vortices. This solution can be regarded in two different 
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ways: as a function with a periodic set of maxima and minima; or, in a more 
physical way, as an ensemble of vortices in plane that interact weakly. Then 
one can easily assume that under the influence of perturbations, the precise 
periodic geometry can be lost and the vortices begin to move rather freely in 
plane and the interactions take the form of inelastic collisions. 

The transformation of ridges having as section the sechyperbolic function, 
into an array of vortices has been known previously and is well documented 
numerically ^Hj- We claim that the actual reason for this evolution, in 
the numerical simulation of the FP equation, is not the attracting nature of 
the monopolar vortex (which actually it is not even an exact solution) but 
the fact that the initial function evolves to the nearest soliton-like solution of 
the closely related equation, the AH equation. Or, this solution is a lattice 
of monopolar vortices. 

We can also explain why the monopolar vortex of the FP equation, ob- 
tained only numerically and studied under a certain numerical imprecision, 
has an amazing stability: it actually represents one only exemplar of vortex 
from the lattice, isolated and extended to 2D infinity by simply ignoring its 
appartenence to the full, periodic solution. This is only numerically accept- 
able, as an approximation, if the distance between neighboring vortices is 
large enough. 

We note, in addition, that the FP equation has also ressemblance with 
the nonlinear differential equation 

( A + A) = sinh (j) 

which describes the stationary states of the Hasegawa-Mima equation. This 
equation, again, is exactly integrable on periodic domains, and its solution 
can be analytically expressed as ratios of Riemann B function. A comparison 
between different functions and the exact solution, represented by p is shown 
in FigEl 

4 Coherent structures and drift wave radia- 
tion 

In previous studies [21], we have developed a theoretical approach to 
describe a coherent structure (the monopolar vortex of the FP equation) 
immersed in a turbulent environment [lOj, |l2], The theoretical in- 

struments are inspired by the quantum theory of fields, in particular the 
calculation of field correlations from the effective action [50], We have 
defined the action functional according to the Martin-Siggia-Rose description 
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Real[Weierstrass P]:blue, lmag[Wierstrass P]:cyan , sech :red, sechigreen 
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Figure 2: Comparison between the sech'^^^ with the Weierstrass function, 
which is the exact solution of the FP equation 



, , and calculated explicitely the generating functional of correlations 

-1/2 



Zj = exp (iS*, 



Js 



2H 



' (2.)-'= ( det '''' 



(3) 



VJs,XJs 



where the determinant is obtained from the product of the eigenvalues of the 
operator 



det 



Ha. 



(4) 



Then the final expression has been derived in the form 

(27r)' 



Zj = exp (iSj) Yl 



det 



<fi.Js,XJs 



-1/2 



(5) 



r /5/2 1 


1/4 


a/2 


1/2 


sinh {(3/2) _ 




sin (a/2) _ 





= const exp {iSj 
and the correlation results from functional derivatives. 



5^Zj 



i6J{y2)i6J (yi) 
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This theoretical approach may be very useful in the study of the late stages 
of the perturbed flow, where the potential consists of both coherent vortices 
and drift wave radiation. 



5 The analytical solution of the FP eq. 



In a previous work jzS] we have determined an exact solution of the FP 
equation. The method consisted in looking to the trajectories of the one- 
dimensional solution singularities, in the complex plane of the spatial vari- 
able. The form of the solution has been assumed as 

Af oo ^ 

and differential equations (plus a constraint) have been derived for the po- 
sitions of the complex poles yn{x). An exact solution to the Petviashvilli 
equation with constant coefficients a and (3 is 

(I)ix,y) = ^ + sp(^iay + ibx + uj\g2 = ^^^ (7) 

with the condition 

a' + b' = f (8) 



where p is the doubly periodic elliptic Weierstrass function. 



6 Comparison with experiments 

Experimental measurements of the characteristics of the zonal flow have been 
performed on Doublet III-D tokamak (|2Z|)- In Ohmic and L-mode plasmas 
it has been found a perturbed potential (prms > 10 and a flow shear ujexb ^ 
2 X 10^ s'^. The value for the radial wavelength is in the range krps G [0.1, 0.6] 
which means G (15. ..30) ps- Due to the different sensitivity of the solution 
(ps to the parameters, this set is sufficiently restrictive to determine its form. 
We have to take v^/u very close to unity and ~ —1500 which gives: 
^rms > i7V, \r — 17. Aps, uexb ~ 2.2 X 10^ s~^. The relative amplitude 
of the perturbation results ~ 4%. An important experimental result is the 
radial spectrum of the perturbation. We have calculated S [kr) from the 
Fourier transform of the correlation of (ps [x). The result is very close to 
Fig.3 of Ref . fT7\ . with two symmetric peaks at /cor ~ 4cm~^ The same 
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Figure 3: Exact solution of the FP equation. 
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Figure 4: The realation between the amphtude of the potential of the flow 
and the width of the periodicity layer, parametrised by the ratio v^/u. 
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The spectrum of the stationary potential perturbation 

0.1 I , , , , , , , 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 O.B 1 

Figure 5: The spectrum obtained from the Fourier transform of the correla- 
tion of the solution of the FP equation. 

sharp decay for \kr\ < k^r is observed, as described in Ref.j2ZI. We also 
note that the analytic spectrum is structured with very narrow dips on its 
range k^. G [1,10] cm^^ which, according to Ref.j2Z] cannot be resolved in 
the measurements. Again, a small fluctuation of the parameters generate an 
average spectrum very close to Fig. 3 of the mentioned work. . 

7 Comparisons with results from numerical 
simulations 

7.1 Numerical simulations at Lausanne 

The conditions are the following: 

= a/3 (m) 

Bt = 2.5 (T) 

Qi = 120 (MHz) 

Ps = 1.8 X 10"^ (m) 

Te = T, = 5 (KeV) 

In these conditions the results obtained were presented in several plots in 
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Ref.jSOj. 



7.1.1 Geometry of the zonal flow 

The radial electric field is shown in Figure 4 of this reference. The values are 
in the range 

Er ~ (-0.02- ■■ + 0.02) X 1.2 X 10^ (l^/m) 
~ (-24 • ■ ■ + 24) X 10^ (V/m) 

Another characteristic of the flow obtained numerically in this reference 
is the periodicity width, 6x. This may be estimated on the basis of the Fig.4. 

(0.06 ■■■0.1) a 

and a in this calculation has been taken 

a = 135pL 

which gives 

6x ~ (8.1^^^13.5)pl 
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Using our analytical solution we can compute the radial electric field on 
the basis of the analytical solution, 0^ (x) . This gives 

^phys ^ 



Tel ( dq.. 



From this we can calculate the flow velocity 

J^phys 



7.1.2 The shearing rate 

In the reference jSO] , the shearing rate is computed as function of the radial 
coordinate s = r/a and represented in Fig.5 and in Fig. 11. 

d flEr 



uJexB = r— . „ 
dr \r B 

The shearing rate is oscillatory as expected for the quasi-periodic geometry 
of the flow, and a typical value at maximum is 

ujexb ~ 2 X IQ-^rii 

We can calculate this quantity 

„,phys J^jPhys 

jphys ^ ^ "^y 

ExB ^phys drP^y^ 

In the slab geometry assumed in our model, 

flo,phys 

phys _ "-"y ( -1\ 

^ExB - ^^phys \^ I 

which can be normaised to fij. 

The results and the particular choice of parameters compatible with the 
simulation Lausannel are in figures. 
Physical parameters 

rhos (m) = 0.0018 

cs (m/sec) = 0.21600E+06 

Omega_i (sec**(-l)) = 0.12000E+09 

Lji (m) = 243.0 

L_T (m) = 0.073 



16 



d(Ln)/dx ( ) =6.5 
Vstar (m/sec) = 0.160E+01 

Parameters in the equations (normalised) 

ekn (rhos/Ln) = 0.0000074 

ekt (rhos/LT) = 0.0246914 

eknp (rhos**2*d/dx(l/Ln)) = 0.0000000 

V* (v*/Omegai*rhos) = 0.741E-05 

u (u/Omegai*rhos) = O.lOOE-04 

cs (cs/Omegai*rhos) = O.lOOE+01 

Parameters in the Flierl-Petviashvili equation 
alpha = 0.259259 

beta = 1.783265 

s = 0.005000 

g2 = 2536.4 
g3 = -2704.2 

Period omega = 0.265228 

Period omega_prim = 0.257673 
Orientation of the stationary flow pattern 

acoef = 0.0000 

bcoef = 0.0385 
Width on x of the periodicity layer 

deltajx (in units of rhos) = 13.3685 

Space domain : 

y (poloidal) ymin = 0.0000 ymax = 100.0000 
X (radial) xmin = -2.0000 xmax = 106.9479 

The minimum of the solution PHIM = 0.727002E-01 

the maximum of the solution PHIX = 0.107068 

The difference of min/max amplitude, si = 0.343680E-01 
The average amplitude, sO = 0.898842E-01 

eddy turn-over time tau (sec) = 0.433340E-04 
Electric field on x, in V/m, minimum : -23237.9 

maximum : 23238.0 
Shear rate normalized at OMCI, min: -0.272078E-01 

max: 0.136097E-01 
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The spectrum of the stationary potential perturbation 
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Figure 8: The flow shear rate and the spectrum S{kx) 

8 Linear stability of the stationary periodic 
solution 



8.1 Linear dispersion relation 

We start from the equation derived before, with the presence of a temperature 
gradient 

|(l-Vi)^=|[(-Vi + 4,^),]-,| 
To study the stabihty we wiU perturb the stationary solution by a smaU 
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function e (x, t) 

(t)^ (t)s ix,y) +e{x,y,t) 

where (ps {x, y) is the periodic flow solution (expressed in terms of the Weier- 
strass function). We obtain 

d d ( 1 

— (1 - A) £ = — -A£ + 477'£ - - -e' 



dt dy 

after taking into account the fact that (f)s is a solution of the stationary 
equation. 

This equation must be solved starting from an initial condition e (x, |/, t = 0). 
If we neglect the quadratic term we have the equation 

d , . s d ^ , ^ , \ de dd>s 

_(l_A).+ -A.-(V-^.)-+e^ = (9) 

If we consider the stationary periodic solution as being aligned along the 
poloidal direction, y, then the function (pg has no y dependence and the last 
term can be discarded 

g^(l-A).+ g^A.-Vg^+^.^=0 (10) 

In {kyjLu) space 

1 /"~ 

e{x,y,t) = 2 / dkyduj exp {—iujt + ikyy) e {x, ky,uj) 

(27r) J-oo 

the equation is 



1 + kl- 



^ {-ico) e + iky (^-k^ + e - iky (477^ - 0,) £ = (11) 



y Q^2 

Some rearrangement yields the equation 



dx^ iuj + ik 
Let us introduce new notations 



1 



1 + ^ 



^ 7 7 2x2 

oi = ujky - ky-4r) 

At,, 
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and the form of the equation is 

+ (giOi + qi(ps)£ = 
The stationary periodic solution is expressed as 

(f)s (x, y) =(f)o + sp (x) 

where 0o is 

and s is a parameter related to the amplitude of the potential perturbation. 
The equation becomes 

0^ + [qiai + qi(l)o + qisp (x)] e = (12) 

This is a Calogero-Moser problem but for a single particle and the analysis 
of this Schrodinger-type problem is, to our knowledge, not available (see 
arxiv.org/hep-th/9903002). In order to obtain an estimation for physical 
cases we will retain the periodic character of the Weierstrass function by the 
replacement 

p (x) ^ So + si COS {px) 



where 



p = (Sx) 

■50 = ^ [P (x) Imax + P (x) |min] 
Si = p{x) Imax - p{x) |mm 



The parameter p is the inverse of the periodicity length accross the layer 
in the x direction. The latter is obtained as 6x from the imaginary-axis 
semi-period of the Weierstrass function and the physical parameters. The 
parameters sq and si represent respectively the average amplitude of the 
function p (x) and respectively the variation of its amplitude on a period. 
This approximation transform the equation (fT^ in the Mathieu equation 

+ [a2 + &2 cos (px)] e = 



dx"^ 



with the notations 



02 = gifli + gi0o + gisso 
62 = gissi 
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We then change the space variable 

, _ px X 



X — > X 



2 2 {5x 

and arrive at the final form of the equation 



dx"^ 

with the new notations 

\2 



+ [a - 2gcos(2a;)]£ = (13) 



a = 4 {8x) qi (ai + 0o + sso) 
q = —2 {Sx)'^ qiSSi 



In detailed form 



a = 4 (5a;)^ ^^^^^ (^-^ - '^K - - Arf + 00 + ss^ (14) 
? = -2(5a;)'ssi— i— (15) 

1 + UJ/ky 

In order £ (x) to be a periodic function on x it must be fulfilled one of a 
discrete set of dispersion relations. We take the simplest approximative form 

% (9) « (16) 



This leads to the following relation between the parameters 

1 

where 

P = -4r]^ + (f)o + sso (18) 
Q = liSxfs'sj 

Eliminating the numitor and grouping the terms 

2 

-{ky + cuf + ^{P-l)-^ + P + Q = Q 

Ky Ky 

We note 

Oi) 

Hy 
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and have the equation 



or 



w'^ + 2tiw + t2 = 



-2kl + P-1 




2ti 



-kl + P + Q 



8.2 Approximations based on evaluations of terms 

Since the spatial variables are measured in ps we can already say that 
wavelength of the y-perturbations will be of the order of several units 

A± > (few units) pg 

We will consider that this is a sufficient indication to take 




Then we can simplify the expressions 




P- 1 



2 




1 - 



P + Q + 1 

1-2 



and the discriminant of the second order equation is 




P- 1 



) 



2 



p + g + 1 

h.2 



A = 



1 - 



2kl 
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The solutions are 



w 



- 1- 



P-1 
'2kf 



P + 1 ± 2kyQ^/'^ 



2kl 



The dispersion relation is at this moment 



Let 



The equation is 



or 



u = 



UJU = — 1 + 



P + 1 



u 



u^ + 2u 



-co ± Q'/^ 
(P + l)/2, 



This leads to 



1 CO ± 



ky 

We can evaluate 
P = 



P + 1 



It 



1 + 



p + 1 



2 (P + 1) 



1/2' 



-477^ + 00 + sso 

-4r7 + -00 + -Sp [X) Imax + 7700 + T^Sp {x) 



+ 



s|mm I f^Vslmax 



where we can have an estimation for the average level of fluctuation 

~ 

no 

In the numerical experiments the field that clearly evolved to breaking 
discrete monopolar vortices had an amplitude of the order of 



0.05. ..0.1 



23 



We chose a case with 

477^ = 1 - — - 1 - 0.85 = 0.15 
u 

The two terms are of similar magnitude so we can expect P to be positive 
or negative. In this case 

P = -0.15 + 0.1 = -0.05 



We take data from one of the characteristic runs 



505 



and 
then 



0s I max 0s | min 

~ 0.088 - 0.019 
= 0.0 69 

Sx ~ 14.17 



Q ~ (14.17)^ X (0.069)^ 



Since 



0.956 



P+ 1 ~ 1 

Q < 1 



We note that Q can be smaller for narrower strips of the stationary solution. 
We can do the following approximations 



2CP+1) 



Q 



1/2 



The dispersion relation is 



1_ _ ujtVQ 

P+1 



It 



VQ 

V2 



1/2 



Q 
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The frequency ui is in the range given by the inverse of the typical time of the 
change of the stationary solution, after being perturbed. The units which 
connect an interval of time in physical and respectively adimensional form 
are ^ 



phys 



with 

~ 0.001 (m) 



ulphys ~ 'f* ~ 0.3 X 10 (m/s) 

0.001 
300 ' 
~ 3 X 10-^ At 



In the numerical simulations the first changes took place after about 

AtUys ~ 4000^7^ {s) ~ Q33^°°°^Q9 = 11940 x lO"'' = 11.94 x 10"^ 
which means 

At\phys 11.94 x 10-6 
3 X 10-6 3 X 10-6 
Then the frequency associated with the change is 

c^^ 1/4 = 0.25 

If we can neglect uj compared with ^/Q we have 



ky P+1 \ ' ^ 

V2 



or, chosing the — sign 



ky 



which has the character of an eigenvalue selection, relating the spatial exten- 
sion of the vortices into which the periodic solution decays to the amplitude 
and radial width of periodicity. 
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Taking from above 
we obtain 



{6x) (50) ~ 1 



X± ~ 27r X 0.41 
~ 2.6 

Our numerical simulations only allow to identify the process of generation 
of a chain of vortices, but we cannot pursue too far the simulations due to 
accumulation of numerical errors. It apparently supports a dimension of 
about 6. 

The stability depends also on the type of perturbation, not only simply of 
the amplitude. A monopolar solution is very stable. The HM term renders 
the equation fragile to decay into lattice of vortices. 



9 Numerical studies 

9.1 Equation with the scalar nonUnearity (time depen- 
dent FUerl-Petviashvih equation) 

9.1.1 Introduction 

The equations ()A.44jl . ()A.47|1 . ()A.46jl . ()A.45jl . help us to define the space and 
time domain of physical parameters for the numerical simulation. 

For the numerical treatment we start form the time- dependent equation 

|(i-vi),^|-(-vi + V),-,| 

and notice that it may be preferable to use two fonctions 

(1-Vl)0 

and 0. The equation is transformed into a system of two differential equations 
for two unknown functions ip and (p 

-vi0 = -(f) + ip 

dt dy ' dy dy 

or 

-vi0 + = i) 



dip _ ^"0 _^ / 1 _^ 4 2\ 

dt dy ^ ' dy dy 
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This version is one of the form taht will be used in the numerical integra- 
tion of the equation. 

9.1.2 Initial monopolar perturbation 

We start by initializing the potential with a monopolar vortex localised in 
one of the layers of periodicity. The amplitude of the potential at this run 
was 0s|max = 0.03 and of the perturbation is 0perf|max = 5 x 10^^ (normalised 
as |e| 0/Te). The evolution is done with a time step of 10f2j . 

We note a very stable evolution of the total flow, consisting of a slow 
[v^^/u < 1) displacement of the monopole along the layer. In the first part 
of the evolution, the profile of the monopolar vortex is adapted to the local 
speed of the flow, but afterwards the evolution is amazingly stable. The 
duration of this period of stability is greater than 12 x lO'^fij the limit where 
the run was stopped. 

This conclusion must be useful in relation with the long stability of the 
Red Spot of the Jupiter atmosphere. There, the configuration is precisely 
the same as in this run: a localised vortex in a structure of zonal flow, as 
shown by the observations of Voyager 1 and 2 (see Busse [2^1 )• 



(|) perturbed, lime = (|i perturbed, time = 29 n s 




10 20 30 40 50 60 70 10 20 30 40 50 60 70 

y/pg (pololdal) y/p^ (poloidal) 



Figure 9: The initial perturbation is monopolar and located symmetrically 
in a flow layer. The initial {t = 0) and final {t = 29/is) potentials are shown. 

9.1.3 Initial dipolar perturbation 

The initial field is (ps perturbed by a dipolar vortex placed with maxima 
approximately on a line of maximum. The amplitude of 0^ is 0.03 and that 
of the dipole is 3 x 10~^. The evolution in this case is very stable for all the 
interval of the run (29 /is, corresponding to 12 x lO^Qi). We notice the weak 
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change of the shape of the dipole due to the interaction with the background 
flow, in the first part of the run. 




Figure 10: A dipolar perturbation is initialised. The initial {t = 0) and final 
{t = 29/is) potentials are shown. 

In another run, the initial field (ps of amplitude 0.03 has been perturbed 
by a dipolar vortex with the same localisation as before but with amplitude 
comparable to the background field, 0.03. The run evolves without destruc- 
tion of the configuration. 



(|i perturbed, dipolar, time = i-i s 



i\> perturbed, dipolar, time = 29 i-l s 



10 20 30 40 50 60 70 

y/p (poloidal) 



Figure 11: An initial stationary flow of higher amplitude is perturbed with 
a dipolar vortex of comparable amplitude. Only the scalar nonlinear ity is 
retained. 



From the numerical experiments with the dipolar perturbations it seems 
that low amplitude potential (pg has weak interaction with the perturbation 
and major changes of the flow may probably arise later than several thousends 
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9.2 Equation with scalar and polarisation drift nonlin- 
earities (time dependent Flierl-Petviashvili equa- 
tion with Hasegawa-Mima term) 

9.2.1 Initial monopolar perturbation 

In a series of run we take a higher initial background flow <^s; with amplitude 
of 0.06 and perturb it with a monopolar structure strongly elongated in 
the radial direction. The amplitude of the parturbation is very high, 0.12. 
We note that the evolution leads flrst to strong changes of the perturbation 
structure which is rotated and deformed. But in the same time the flow 
itself tends to be broken into isolated vortices. These runs require higher 
numerical precision so that only short period of time can be afforded at this 
momant. The evolution was stopped at 6 fis. 



i|) perturbed, transversal, time = !.i s (|i perturbed, transversal, time = 6 1.1 s 




111 20 30 40 51) 60 70 10 20 30 40 50 60 70 

y/p^ (poloidal) y/p (poloidal) 



Figure 12: A run with both scalar and polarisation nonlinearities retained. 
The flow is initialised with a monopolar form perturbing the staionary peri- 
odic solution. 

The late stages of these runs are useful for the ilustration of the breaking 
of the background flow into a lattice of vortices. 

A particular importance takes the examination of the evolution of the 
flow when the initial perturbation is monopolar and placed approximately in 
a layer of periodicity. This should be directly compared with the previous 
runs with the pure scalar nonlinearity. 

We notice that the stability is much stranger for this type of geometrically 
"compatible" perturbation. Neither the background flow nor the vortex is 
disturbed for more than 25 /xs. The monopole is only advected by the flow. 
At late stage, the monopole is deformed and destroyed, while the flow does 
not show the previously observed tendency to break into a lattice of vortices. 
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(|) perturbed, transversal, time = 6 |.i s i|] perturbed, transversal, time = 6 1.L s 




10 20 30 40 50 60 70 HI 20 30 40 50 60 70 

y/p^ (poloidal) y/p^ (poloidal) 



Figure 13: The late stage of a run with monopolar pertubation.The breaking 
of the flow into vortices is apparent. 

However, the region of the very high amplitude where the monopole existed 
tends to break into vortices. 



FP&HM (ji perturted, monopole, time = n s FP&HM <ti perturted, monopole, time = 27 ).i s 




10 20 30 40 50 60 70 10 20 30 40 50 60 70 

y/p^ (poloidal) y/p^ (poloidal) 



Figure 14: The late stage of a run with monopolar pertubation of high am- 
plitude. The breaking into vortices affects first the strongest perturbation. 

This shows that the Hasegawa-Mima nonlinearity has the property of 
inducing a tandency of the flow to break up, under strong perturbation, into 
a lattice of monopolar vortices. 

9.2.2 Initial dipolar perturbation 

A dipolar perturbation shows essentially the same evolution as for monopolar 
case. The profile is not changed for rather long period, 20/is but the later 
stages consist of a clear decay into a set of distinct formations with monopolar 
structure. We cannot follow too far the simulation due to accumulation 
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of errors, but it is clear that space scale of the new vortical structures is 
close to those found in preceeding cases. This strengthens the idea that the 
polarisation drift nonlinearity is the cause of this evolution toward lattices 
of vortices. 



FP&HM (t) perturbed, dipolar, time = fi s 



FP&HM [|j perlui-bed, dipolar, time = 28 ,u s 





Figure 15: Initial dipolar perturbation 



9.2.3 Initial oscillatory perturbation 

In order to examine the processes of destruction of the flow and decay into 
isolated vortices, we initialize with an oscillatory perturbation with a period 
much larger than that which seems to be chosen by the system as extension 
of the final monopolar vortices. The flow clearly evolves to the lattice of 
vortices. 



FP&HM (|i perturbed, space periodic, time = i^l s 



FP&Hlv! (|> perturbed, space periodic, time = 9.5 u s 





y/p (poioidal) 



y/p (poioidal) 



Figure 16: Perturbation initialised as a spacial oscillation of large wavelength. 

For a higher amplitude of the initial perturbation the evolution is quali- 
tatively similar. 
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Figure 17: Initial oscillatory perturbation at higher amplitude. 



In a series of runs, we initialise with a perturbation consisting of oscilla- 
tions with higher spatial wavenumber, i.e. slightly closer to the final vortices. 
The result shows clearly that the initial conditions do not influence signifi- 
cantly the space extension of the final monopolar vortices, which means that 
the system chooses this dimension according to a dispersion relation which 
has the character of an space-eigenvale problem. This is also suggested by 
the fact thet the final stages appears to be similar for rather different time 
duration which is needed by the system to reach the stage of breaking into 
vortices. 



FP&HM perturbed, many space periods, time = m s FP&HM i|] perturbed, many space periods, time = 7.7 yi s 




Q: 10; ^ 30 40 50 60 70 10 20 30 40 50 60 70 

y/Pg (poioidai) y/p^ (poioidal) 



Figure 18: Initial perturbation taken as a space oscillation of relatively small 
wavelength. 
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9.3 Conclusions from the numerical studies 

Extensive numerical simulations (which however must be considered prelim- 
inary) have been done for studying the stability of 4>s- In a series of runs 
the time dependent equations with pure scalar nonlinearity has been used 
with (f)g perturbed by (f)p, taken first as a monopolar vortex . For small am- 
plitudes {e.g. (j)s ~ 0.03 and 0p ~ 0.005) the total flow is stable for more 
than 12 x 10"^i7j. The monopole is reshaped at early stages for equality of 
the tangential flows and then is stably advected {v^/u < 1). This result may 
be relevant for the stabihty of the Red Spot, the long-lived vortex embedded 
in the zonal flows of the atmosphere of Jupiter, a geometry very similar to 
ours. A small amplitude dipolar initial perturbation is also stable for similar 
durations. 

The structural stability of 4>s is investigated using the time dependent 
equation with both scalar and polarisation drift nonlinearities retaind. The 
flow is stable on only shorter duration and a new effect arises: at late stages, 
the flow tends to break up into a set of monopolar vortices. Taking (f)p of 
high amplitude (0^ ~ 0.06 , (f)p ~ 0.12) monopolar and strongly elongated on 
r there is a clear tendency to rotate it such as to align with the background 
flow. At 6 /IS the process of generation of vortices is pronounced. Similar 
conclusions can be drawn from runs with 0p taken as dipolar and as spatial 
periodic structures. Examination of many runs leads to the conclusion that 
the spatial extension of the generated vortices depends only weakly on the 
initial conditions and of the time duration of quasi-stability and is A ~ 6ps . 
Higher numerical precision is required to study their evolution. 
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10 Appendix 



10.1 Derivation of the equation with time dependence 

The aim of this Appendix is to review the derivation of the various forms of 
the equations used in this work. It will be reviewed the stationary and the 
time dependent equations, with scalar and polarisation drift nonlinearities, 
etc. The main practical purpose is to compare the units used and to facilitate 
the present numerical computation and the future possible changes. For this 
reason some calculations are trivially simple and may be skiped. 

Consider the equations for the ITG model in two-dimensions with adia- 
batic electrons: 



+ V-(vini) = 



dt 

^ + (v,-V)v, = A(_V0) + — v.xB 
ot rrii rrii 

We assume the quasineutrality 

Tli ~ ?T/g 

and the Boltzmann distribution of the electrons along the magnetic field line 

|e| (f)^ 



He — no exp 



e 



In general the electron temperature can be a function of the radial variable 

The velocity of the ion fluid is perpendicular on the magnetic field and is 
composed of the diamagnetic, electric and polarization drift terms 

= ^dia,i + + Vpol^i 

Ti 1 drii^ 
\e\ B Hi dr ^ 
V0 X n 



B 

1 f d 



\dt 



7^ + K ■ v^) 
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The diamagnetic velocity will be neglected. Introducing this velocity into the 
continuity equation, one obtains an equation for the electrostatic potential 

Before writting this equation we introduce new dimensional units for the 
variables. 



>phys 



(^.phys „phys\ 
(A.2) 
Ps ps ) 

^Phys ^ ^/ _ ^phys^, ^^ 3^ 

The new variables {t, x, y) and the function are non-dimensional. In the 
following the primes are not written. With these variables the equation 
obtained is 

I (1 - vi) ^ 

- (-V_L0 X n) ■ V* 
+ (-V_L</) X n) ■ vt0 
+ [(-VX0 X n) ■ Vx] (-Vi0) 
= 

where 

V, = -Vxlnno- VxlnTe (A.4) 
vr = -V±lnTe 

(This is Eq.(8) from the paper [Zj). 
NOTE 

Before continuing we compare this equation with the equation of paper 
|29j . Eq.(16). Here taking still the units to be physical, the form of the latter 
equation is (Eq.(12) from that paper) 

d \e\(h d 1 _o , 



dt Te dtBVti 

H ^5 Vxlnwo 

JD 

— V_L0 X n ^ |e| 



(A.5) 



B T, 



e 







^ [(-Vx0 X fi) • V J (-V^ 
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The term containing the gradient of the equihbrium density comes from the 
continuity equation, as convection of the equihbrium density by the fluctu- 
ating E X B velocity. The term containing the gradient of the equilibrium 
temperature comes from the continuity equation by the divergence of the 
fluctuating flow of particles, where the adiabaticity has been assumed. 



The equation of jSHl remains expressed in these terms. However, this is 
strange in relation with the derivation made by Lakhin et al, where the 
temperature term is canceled. See below a discussion of the Eq.(8) from 
Laedke and Spatschek 1986. 
For the second term we have 



n \e\ (j) 



no Te (x) 




BVti \e\ ^ Te 




This will become (with its sign) 




(A.6) 



in the new variables Eqs. ()A.Hl - (|A.3jl . 
The third term is 





e| (j) 



X n ■ V_L In no 




This will become 



i-Ps'^L^) X n) ■ (-psVxlnno) 



in the new variables 



a(-Vl(/)' X n) ■ (- VI In no) 



(A.7) 
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The fourth term is 



X n 



1 



B 



rrii 



-V 



\e\cj> 



X n 



\e\l 
e\4> 



X n ■ V J 



^iPs -V±^ X n 



Te 



X n 



X n 



7 |e|</> 

-'■ e 



|e|0VxTe 



Te Te 



T. 



V±Te 

T. 



This term will become 



(-V'^c/)' X n) 4>'-{-VjnTe) 



(A.8) 



after introducing the new units. 

The last term (with the polarization nonlinearity) is in physical units 



1 



[(-Vx0 X n) • VJ (-Vi0) 



This is converted to non-dimensional variables 



1 



Ps |e| 



Hi 

Te 



X n 



1 



pI |e| 



Te 



Collecting the physical coefficient we have 



Te 



1 



1 /Te 



\ mi J 



Pi 



Then, in the new variables, this term becomes 

[(-V^0' X n) • VJ (-V20') 



(A.9) 
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Then the Eqs. flA.5|) with the new form of its terms ()A.6|) . ()A.7|) 
and ()A.9|) becomes 

-n, (-V10' X fi) ■ (-Vllnno) 
+VLi X fi) 0' ■ (-V1 InTe) 

+a [(-Vlf X n) • VY] (-V^f) 
= 

Introducing the time unit and ehminating the primes 

|(l-Vi)0 (A.11) 
-(-V±0 X n) ■ (-V±lnno) + (-V±(/) x n) ■ (-V±lnTe)0 

+ [(-v^0 X fi) ■ v^] (-vi0) 

= 

However this is not the equation (16) of Laedke and Spatschek 1988. This 
is discussed in the following. 
END OF THE NOTE 

We return to the physical equation, Eq. ()A.5j) or Eq.(12) from [221 (1988) 

d \e\(l) d 1 _o , 



dt % dtBQi 

H 7^ Vxlnno 

— V_L(/) X n ^ |e| 



(A.12) 



= 

We make a change of variables 

7] = y — ut 

which means 

d d d 
dt dt' dr] 
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The equation becomes, after removing the primes from t 

d \e\(j) 9 1 , d leld) 9 1 , 
-Vi0 - M— ^ + M— — — Vi0 



dt' Te dt' BVli dfj drj BVti 

H 5 V±lnno 

is 

, -Vx(/) X n _^ |e|(/) 



e 



= 

The first term is 



(9 |e| (9 1 |e| -B 



The third term is 



dt % dtB mi Te 

d 1 VtiVti 

1 a 1 

BQidi^/ 



9 lei 1 9 1 

—u— — — — = — „^ u- 



(A.13) 



dr] Te B^i dripj 
the second and fourth terms remains unchanged 



Since the product 1/ (JljS) arises systematically, we will extract it from the 
third term 

-V±0 X n I Q.i d<j) 

V±lnno = —BFr^—>^n^ 

B B\li u or] 

where we define 

Kn^x = V^lnrio 

with Kn measured in m'^. 

The term with the gradient of temperature 

■ V I — — = e • -—V I T, 



B Te B ' V T? 

after introducing the definition 

Krex = V± In 
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we obtain 



-V±0xn ^ \e\(f) 



B 
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kt 1 |e| B dcf) 
111 dcj) 



The term with the polarisation nonlinearity is 



[(-Vx0 X n) • VJ (-V 



The form of the equation is 

1 d 1 



d 1 



Bn, dt pi ' dt BQi 



T(i) 



Vi0 



19 1, 1 90 

1*7^ — ^T^ti 



BQ,i dr] p1 BQ,i u dr) 
d 1 _o 



+u 



dr] BQ, 



111 d(t) 

"'mi-piB^^^ 



i- [(-V^0 X n) • V J (-V10) 



= 



Multiplying by fli 



d_ 
di 

—u 



1 ^ ^i'^n 



d_ 
dr] 



+ - 



1 Kt 



B 



plB^^dT] 

[(-Vx0 X S) ■ V J (-v: 

= 



(A. 14) 



NOTE ON THE COMPARISON WITH Su and Horton 
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We can compare this equation with the corresponding one, Eq. (45), from 
Ref. [?]. This equation is 



72 \ 9^ 



T x) 



+ [VdO + VdO^ - f^T^) ^ 

- [(-Vxv? X n) ■ VJ {-Vlif) 



where 



T(x) 



To 

Ln |e| $ 

PsO 



This is the small parameter of the expansion. 

Vd = l + v'^qX 

I dvd 

^dO = PsO-^ ~ £nVdO 

END OF THE NOTE 

In Ref.pnj the potential is redefined as 



^ B 



and the equation becomes 



d_ 
di 



U 



d_ 

dr] 



4 + ^1^ 



pt 



u 



drj p1 B"^ dr] 

+ [(- Vx^ X n) ■ V J (- Vl^) 



In Ref.j22j| the following operator is introduced 



P 



1 



2 -J- 



P 
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and the equation is expressed as 
d 

PqI^' = [V± - ux) X n] • {P^ + ^ilnno) 

In this way the term containing the temperature gradient Kt is hidden in the 
application of the operator V_l on 1/ pi contained in P. 



10.2 The stationary form of the equation 

The most general form of the solution at stationarity is 



'dt 



= 



Pijj + Qj In no — 0,i In 



no X 

u 



In the natural logarithm we expand the density at x and then parform 
an expansion 

In (1 + PS £ — — 



Q,i In 



no \ X 

u 



Q,i In 

fli In ^no {x) 
Qjln [no (x)] 
+a In 



, . ip dno (x) 1 ip"^ d^no {x) 
^ u dx 2 dx^ 



^ il) 1 dno {x) 1 d^no (x) 

uno dx 2 no dx^ 

^ ijj 1 dno (x) ^lijj^ 1 dj^no {x) 
uno dx 2 V? no dx"^ 



In 



no\x 

u 



= Vti In [no {x)] 

1 dno{x) 
u no dx 

The last two terms are 

1 1 d'^no (x) 

2 no dx^ 



1 ■ijj'^ 1 d^no lifj"^ I f dn, 

2 no dx"^ 



2 nl 



11/ dno 



2no \ dx 



l_d_ 
2dx 

2 " 



da; 
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The final form is 



Qi In 



no[ X 

u 



= ili In [no (x)] 

In tflis expression tlie units are pliysical: 

u {m/s) 



Tlie equation is tlien 

Pip + VLi In rtQ = Qi In 



no a; 



or 



iljln [no (x)] + ilj Kn + TT^^n 

/ u 2 



(A.15) 



10.3 The equation with time dependence and polari- 
sation nonhnearity retained 

In order to write an equation for the potential with time dependence, we 
return to Eq. ()A.14|) . We multiply by B 



d_ 

dt 



Pi 



-u 



d_ 

dr] 
d 



1 



u 



1 Kt ,(- 



+ 



pI B drj 



+^ [(-Vx0 X fi) . V J (-Vi 




(A.16) 
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Then we derivate Eq. ()A.15|) with respect to rj, multiply by u and return from 
tp to 4> = ijjB which will be noted (pg (stationary solution) 

+ ^^^n^<Ps-^^^ = (A.17) 
or] uB or] 2 



1 ^ . ^2 d 
or] or] 



Now we add this equation from Eq. ()A.16|) and combine the two functions 
and in order to obtain an equation for their difference. 
We first examine the last term 



{[-Vx(0-0s)xn]-V^}[-V^ 
= [(-V^0 X n) • V J (-Vi0) - 

- [(-Vx0s X n) ■ Vx] (-Vl0) 

- [(-Vx0 X n) ■ Vx] (-Vi0,) 
+ [(-Vx0s X n) • Vx] (-Vl0,) 

The last term is zero. Expressing the term we need (the first one) 

[(-Vx0 X n) ■ V J (-Vl0) 
= {[-Vx(0-0s) xn]- Vx} [-V^ 
+ [(-Vx0. X n) ■ Vx] (-Vl0) 
+ [(-Vx0 X n) ■ Vx] (-Vl0,) 

After adding the equations we obtain 



d_ 

dt 



Pi 



d_ 

df] 
d 

+ ^ [-^iK^n (0 - 0s)] 



+ 



1 Kt 0^ 

+^{[-V^ (0 



1 



a<0r 



1 



m5 2 

-0,)xii]-vx}[-vi (0-0.)] 



+ ;g [(-V^0, X n) ■ Vx] (-Vi 




- [(-Vx0 X n) ■ Vx] (-Vi0,) 



We define the new potential representing the difference between the time 
dependent solution and the stationary solution 



99 = - 05 
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Then the last two terms of the above equation will become 

[(-V±0s X n) ■ Vx] (-Vi0, - Vi(^) 
= [(-Vx03 X n) • VJ (-Vi(^) 

and 

[(-Vx0. X n) ■ + (-Vxv^ X n) • V J (-Vl 
= [(-Vx(^ X n) • V J (-Vi0,) 

We can use in the last term the equation defining the stationary solution, (pg- 



1 ^ fijK^ 

p2 u 



u^B 2 



a(t)s - I3(t)', 



with notations 



a = ^ + 



/3 



Then the equation is 



d (\ 



dt 



-u 



5^ VPs 

d_ 



d_ 
dr] 
1 



{(f)s + <p)- 



r0s 



2p25 ' 2uB 
+^ [{-^±V X n) • VJ {-Vlifi) 



B 

[(-Vx0. X S) ■ Vx] (-Vl(^) 

+^ [(-Vx(^ X n) • VJ (-Vic/),) 
= 

In this equation, the quantitities are in physical units 

{x,r)) m 
t s 

(p V 

u m/ s 



(A.18) 
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We change to the units 



Also 

lei 



After some arrangement, we have 



a^(i-v2),' (A.19) 



7 (i - v!) ^' 



d 

I 

yPhys Q 

fliPs dr] 

B 



B\e \ Ps^i I 2 dr]' 

P^^^ ,jphys\ 1 A r^M2 



/2 J,/ 







In this form, the prime in 0^ means that it will be calculated for (a;', ry'), i.e. 
with normalised variables. 

We introduce new units for the speed of the reference system, u 



yPhys 



Also 



u=-— (A.20) 

\liPs 



ri'n — Psf^n 
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physrp ^phys 2 

txirp J. g tvrp C-^ 



Then we introduce 



phys 



phys 

Kt = PsI^T 

fhr^^_ Iphys ^ ]i„2 /phys 

\e\ B uP y u 
where u is the adimensional speed, Eq. ()A.20j) . It is then useful to introduce 

k' — r?K'P^y^ 

and the coefficient in the fifth hne of Eq. ()A.19jl is, in terms of adimensional 
quantities 

PsTe 1 .Jphys _ ^ 



Finally we note that 



T 



Bpl \e\ ni 

The same adimensionalization is done for a and f3. The equation is multiplied 
by p1 and (ps is normalized to 0^. 

\pI up^y' J up^y / {VtiPs) 



1 + ^ 

u 



T n.^'phys rp 

\ey 2{uP^y'fB\ey 



2 VA^. n ; «n2(«P^ys)2/(fi2) 

1 (pyp^y^) 



2 



1 K,'^ 



2v? 

The equation in adimensional variables (and removing the primes from 
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the variables and operators) becomes 

d_ 

dt 



4(1-Vl)^ (A.21) 



d 
or] 

u 2dr]^' 
+ [(-Vx(^ X fi) ■ Vx] (-vi^) 
+ [(-Vx0, X n) ■ V±] (-Vi(/^) 
+ [(-Vx^ X n) ■ Vx] (-Vi0,) 
= 

In contrast to Eg. ()A.18j) this equation is only expressed in terms of adimen- 
sional variables. 

Let us introduce a new function 

^ = (1 - Vi) (A.22) 



Then 



^ _ «^ _ (A.23) 
dt dr] dr] 

+ [(-Vx<^ X fi) ■ Vx] - <^) 
+ [{-V±(ps X fi) ■ V^] - (/?) 
+ [(-Vxv^ X n) ■ V J (-Vl 




The last three nonlinear terms are 



,-V...a).(V.,-V..)^-|g.g|| (A.4) 

;-V.,.xfi),V.,-V.., ^ -f f I (A.25) 

c^v? d(f)s dip 



drj dx dx drj 
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(A.26) 



drj dx 

The final form of tlie equation for is 



s ^ dip d(t)s 



dx dr] 



dip dip d(f 
dt dfj ^ df) 

dip dip dip dip 

dr] dx dx drj 



d(ps dip d(ps dip 



+ (-a + 2l3(ps 




dip d^s ^ dip d^s 



dr] dx dx dr] 



(A.27) 



dr] dx dx dr] 
2dr] ^ u 2dr] ^ drj dx dx dr) 



10.4 The time evolution of a perturbation around the 
stationary poloidal flow 

In this case, cps is the periodic solution expressed in terms of the Weierstrass 
function, and with parameters chosen such that the pattern of the flow is 
parallel with the poloidal direction. Then (ps only depends on x and 

(A.28) 

which considerably simplifies the equation. 

In addition we will take again the fixed (laboratory frame, for this case 

t t' = t 

r] —>■ y = r] + ut 

which gives 

d_ _ d_ 

dr] dy 

d_ _ d__ d_ 

dt' dt dy 



49 



and the equation becomes 



dip dip dip dip dip d(ps dip 
dt dy dx dx dy dx dy 



d(ps dip dip Id 
dx dy ^ dy 2 dy 
dip d(ps 



+ + dydx 




or 



dip ^ dip f dip\ ^dip fdip ^ 



dt dx \ dy J dy \dx dx 



(1 - a + 2l3<Ps) ^ + ^^r. 



(A.29) 



This can be written 



where 



Dip 
~Dt 



+ C[ip\ = ^ 



Dip dip dip 

= h Vt 

Dt dt dx 



dip 



(A.30) 



Vx = - 



dip 
dy 
dip _^ 
dx dx 



(A.31) 



C\ 



(1 - a + 



d(ps 
dx 



-Kn - {1 - a + '2P(ps) + I^T {(ps + f) 



dip 
dy 



The solution of Eq. (fOn|l is 



ip (x, y) =ip {xo, yo) - / drC [ip {xq + tv^, yo + rvy)] 

J to 



(A.33) 
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The Eq. ()A.33|) gives a possible way to solve numerically the equations. 
The final system is 



- Vlv5 + 



dr 



(A.34) 



to 



-K„ - (1 - a + 2p(j)s) + (0, + if) 



dip 
dy 



(A.35) 

y=yo+'i'Vy 



We note that the vector field of Eq. ()A.3H1 does not have zero divergence 



Vi -v 



(A.36) 



10.5 The equation of Petviashvili with time depen- 
dence due to the temperature gradient 

10.5.1 Physical model of the ion drift instability 

This equation has been derived in several papers and used in numerical cal- 
culations (Laedke Spatschek 1986). From the dynamical equations of ion 
density (as above) it is obtained 



(A.37) 



X n 



-V_L0 X n) ■ Vdr0 



[(-Vx0 X fi) ■ V J (-V 



where 



- V_L In no — V_L In Tg 



-V, InT, 



This is Eq.(8) from [?]. It differs of what we have obtained before. 

Later the equation is expressed as, explained in the article as the limit 



oo 



This is the equation (12) from paper of 1986 by Laedke and Spatschek j?]. 
Here the units are: time is measured in space in p^; the potential (f) is 

|e| (t)P^y'/Te. 
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10.5.2 The transformation of the equation 

A referential moving with the velocity u is introduced 



y 

t 



y' = y-Ut 
t' = t 



Then 



or 



1 d 



dy' 2 dy' 



A(i-vi),-4(i-vi), + ..|-„4A,'^^o 



We will drop the prime 
d 



d 



We divide by u 



udt ^ dy 
and replace the time variable by 



(A.39) 



2h 



t^t' = ut 

Here the time in Ea. (l09|l had been already adimensionalized by expressing 
it in units of ■ Since y in the same equation is measured in p^, the velocity 
is measured in units of 

Ps 

This is the case for both u and v^. The equation becomes 



with 



u — v* 



> 



u 



(A.40) 



(A.4i: 



We now renormalise the potential by including in it the factor —Vdr/u. 



(A.42) 
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For this both sides of the equation are multiphed by —v^t/u and the equation 
is 

|(l_Vi).^' = |(-Vi + 4,=)^'-0'|: (A.43) 



The units are: 

e(t>\phys VdT e(t)\phys 



iphys 



0' = 



The velocities are measured as 



VdT e4>\phys 



or (A.44) 



u T, 



e 



U\phys —p^u\phys Or (A.45) 

' - -J— I 

U — ^_iU\phys 

This quantity disappears from the equation, being absorbed into the new 
time variable. The time is measured as 

i\phys ^ ^Iphys ^ ^ ^\phys Or 

t — '^Iphys^i t\phys 



Ps 

The distances are measured in p, 



1 

"^iphystlphys 



Ps 

Naturally, the physical variables are measured in SI 



y' = —y\phys (A.47) 



Vlphys 


m 


1 1 phys 


s 




mj . 


4^\phys 


V 



and at this point the primes are suppressed. 

If we have a result about some time duration At' {i.e. expressed in terms 
of the non-dimensional variables) then in order to recuperate the physical 
duration, we have to do 

u\ 



\phys 



53 



Only for 

u > f * 

the localised stationary solutions are posible. 

The stationary form of the Petviashvili equation is 



|^(-Vi + V) 0-^1 = (A^48) 

A0 = Arf<p - (A.49) 



which is the same as the general form for 
These can be found in 



10.6 The equation derived from Ertel 's theorem 



The derivation of the equation is done on the same basis as the derivation of 
Laedke and Spatschek 1988. However, in a series of papers, Su and Horton 
|H], develop a more physical justification of the scalar nonlinearity model, 
although their primary aim was to investigate the stability of the dipolar 
(Larichev-Reznik) solution of the Hasegawa-Mima equation. 
The equations for the ion instabilities are 



d_ 
di 
dn 
'di 



+ (v ■ V) 
+ V I (nv) 



rrii 



It is derived the Ertel's theorem 



dt 



^» + n ■ (Vx X v) 

n (x) 







(A.50) 



This essentially means that: the sum of the cyclotron frequency and the 
vertical component of the vorticity, divided by the density, is constant along 
the line of Lagrangean evolution with the velocity E x B. 

If one neglects the cyclotron frequency and takes the density constant, 
this theorem is equivalent with the Euler theorem for ideal fluid 



duj 







The ordering is 



1 d v-Vi 



^* dt 



< 1 



54 



The equation becomes 



+ [(-V^^ X n) ■ VJ ( i±%Zk^ = (A.51) 



dt \ n J ' 'V 
It is assumed the Boltzmann density distribution for the ions 

I le| $ 

n — uq {x) exp 



(A.52) 



= no [x) exp 



10.7 Stationary traveling solutions 

This means solutions of the type 

cp = cp{x,y- ut) 

where u is the speed of the system of reference. The stationary equation in 
the moving frame is 



I rr-T I + [(-VxV' X n) • Vj 



1 + enV\^ 



no [X ) exp 



T{x) 



no {x) exp 



T{x) 



This equation is identically verified if the expression in the brackets is a 
general function of 

(p — ux 



I.e. 



1 + EnVlif 



no (x) exp 

One can take this function to be 



T{x) 



F {(f — ux) 



F {ip — ux) = no (x — — ^ 



which gives 



no ix) 



no {x - ^) 



exp 



T{x) 



- 1 



(A.53) 



In the case where the density profile is exponential 

no {x) = exp {-Enx) 
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the equation becomes 



£nVi<y2 = exp 
The dimensional parameters are 



1 1 



T (x) u 



Thaking the logarithm of the Eq. ()A.53|) 

In (l + £„V5_V5) = In no (x) — In no — — j + 
The density profile is expanded 



T x) 



v' 



In no (x) ^ -En [ Vdox + -^x^ + 



gives 



Vd (x) = Vdo + v'^qX 
= ^ + VdoX 



The ordering is 



{AM) 



> dvd , . 

Vdo = Pso-^ ~ £nVdO (A.56) 



In (1 + EnVl^f) ^ Enk^ («, x) ^ + £n^f^ + 



k{u,x) = —— A.57 

I ix) u 



pIo^' ^^-^^KTr~^v',,r^^=En^E<^l (A.58) 
1 p_ Lir, 



The velocity of the moving referential u and of f ^o are of the same amplitude, 
normalized at 1. 

U ~ VdQ = 1 

Keeping only terms of order 



^l^ = kl^ + ^y (A.59) 
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where 



with the notations 



k'^ {u, x) = + ax + 



,2 1 ^do 

kr. — 1 ~ £ 

u 

VdO 2 

a — kt ~ £ 

u 



The gradient of the drift velocity gives rise to a hnear damping 



u 



-x(p 



The temperature gradient generates also a linear term that compensates the 
hnear term of the gradient of the drift velocity. 

10.8 The energy of the solutions 

The density of energy has the expression 



The units are 



Eix,y,t) = - 



(^2 (x) 



T{x) 
L„, lei $ 



Pso To 

and T (x) is adimansional. The distances are measured in psQ. This means 
that 



E{x,y,t) 



PsO 



1 



1 



T{x) V To 



\e\^V |e|$^ ^ 
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